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Abstract. We analyze a sample of over 1000 stars from 
32 fields imaged in the V and I bands with the Wide 
Field and Planetary Camera, on board of the Hubble 
Space Telescope. The fields are located at Galactic lati- 
tudes \b\ > 15° and in various directions on the sky. We 
consider models for the structure of the Galaxy with dif- 
ferent choices for the main parameters governing the shape 
and luminosity function of the thick disk and stellar halo. 
Comparing model predictions with the observed colour- 
magnitude diagram we are able to rule out an increasing 
or flat stellar luminosity function at the low-luminosity 
end. We also rule out large values of the vertical scale 
height of the thick disc, zq, finding it to be in the range 
800 < zq < 1200 pc. As for the local density normal- 
ization, values within the range 4% < no < 8% seem to 
better reproduce the data. Our data essentially rule out 
a flattened stellar halo (c/a ~ 0.5) or models with both 
large local normalization and effective radii. 
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1. Introduction 

Our Milky Way galaxy is known to have a large disk with 
spiral structure, typical of luminous galaxies of Sb type. 
It has also been known for a long time to contain a po- 
pulation of old (t ~ 10 Gyrs) stars distributed in an 
spheroidal component. More recent models of its struc- 
ture have included at least a third component, of inter- 
mediate properties, often called the thick disc (Gilmore 
& Reid 1983). Our understanding of the spatial, kinemat- 
ical and chemical structure of the Galaxy has continu- 
ously improved in more recent years as a result from the 
quest for a more quantitative and detailed picture of how 
the Galaxy (and galaxies in general) formed and evolved. 
However, despite the enlarged and improved stellar sam- 
ples with accurate photometric, parallax and velocity mea- 
surements, several issues concerning the structure, origin 
and evolution of the main components and populations in 
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the Galaxy still remain to be settled (Gilmore et al. 1995, 
Majewski et al. 1996, Norris 1999). 

The very existence of the thick disc as a discrete com- 
ponent, with stars kinematically and chemically distinct 
from those of the thin disk and halo, is subject to contro- 
versy (Norris 1987, Carney et al. 1989, Reid 1998). There 
is little doubt, however, that models with one planar com- 
ponent, with a single density profile, do not provide a suit- 
able description of the structure of the Galaxy, since such 
models do not successfully fit recent star count data (San- 
tiago et al. 1996a, Buser et al. 1999). 

The structure of the thick disc is generally described as 
a double exponential, with horizontal and vertical scales 
in the ranges r ~ 2.5 — 4.0 kpc and 600 pc ~ z ~ 1600 
pc, respectively (Chen 1996, Robin et al. 1996, Buser et 
al. 1999, Norris 1999). The wide range of scale height 
values reflects in part the difficult task of defining and 
separating thick disc stars from those of other structural 
components. The value of zq for the thick disc is anti- 
correlated with the local density of intermediate popu- 
lation stars, tiod, which normalizes the assumed density 
profile: models with larger zq usually have smaller uqu 
values and vice versa. 

The spatial structure, stellar luminosity distribution 
and shape of the outer regions of the spheroid, the stellar 
halo, have also been under closer scrutiny in recent years 
(Wetterer & McGraw 1996, Gould et al. 1998, Stetson et 
al. 1999, Moore et al. 1999, Samurovic et al. 1999). Largely 
due to Hubble Space Telescope (HST) observations, it is 
now known that the stellar luminosity function (LF) both 
of globular clusters and field disk stars decreases at the 
low-luminosity end (Dc Marchi & Paresce 1995a, b; Elson 
et al. 1995; Santiago et al. 1996a,b; Gould et al. 1996,1997, 
Piotto et al. 1996, 1997, Mendez & Guzman 1998). How- 
ever, there is no consensus yet about the shape of the 
stellar LF at the faint end for other components (Dahn 
et al. 1995, Gould et al. 1998). Also, the flatenning of the 
stellar halo and the density profile function that best de- 
scribes it are subject to uncertainty. 

Improved observational data will more efficiently con- 
strain the models if use is made of objective and statis- 
tically sound means of comparing model predictions to 
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Table 1. The new 15 fields. 
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observations. Recent efforts in the direction of efficiently 
modeling observed colour-magnitude diagrams (CMDs) or 
other N-dimcnsional spaces of observables have been de- 
veloped (Saha 1998, Hernandez et al. 1999, Lastennet & 
Valls-Gabaud 1999, Stetson et al. 1999). These methods 
usually try to make use of non-parametric statistics and 
a minimum amount of initial assumptions. 

In this paper we analyze a composite CMD of 32 stel- 
lar fields towards various directions in the Galaxy, all im- 
aged with HST's WFPC2. We apply simple and objective 
statistical tools to compare the observed CMD with those 
generated using structural models of the Galaxy. In Sect. 2 
we discuss the data used, whereas in Sect. 3 we present the 
statistical methods and models. In Sect. 4 we discuss the 
models found to best reproduce the data and the consis- 
tency among the different statistics used. Finally, in Sect. 
5 we present our final conclusions and future perspectives. 

2. The WFPC2 data 

Our WFPC2 fields are part of the Medium Deep Survey 
(MDS) database. The data were extracted from 32 deep 
fields, 17 of which, mostly at high Galactic latitudes, were 
analyzed by Santiago et al. (1996a). The new 15 fields 
are listed in table 1. Their MDS id is given, along with 
their Galactic longitude and latitude (/// and bn, respec- 
tively), total number of exposures and exposure times and 
the number of stars found in each. All fields have been ob- 
served with the F606W and F814W filters and have at 
least 2 exposures per filter. Several of the 15 new fields 
have \bji\ values in the range [15°, 30°], thus containing 
a larger fraction of thick disc stars than those studied by 
Santiago et al. (1996a). 

Sample definition, photometry, instrumental and com- 
pleteness corrections, and calibration to the standard sys- 
tem for the new fields were carried out in a similar fash- 
ion as in Santiago et al. (1996a). In brief, the DAOPHOT 
package within IRAF was used to automatically make an 
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Fig. 1. a) 8m vs Ifsi4 plot for one of our intermediate 
latitude fields, b) \ 2 vs Ifsu plot for the same stars on 
panel a. Indicated in both panels are the cuts done in order 
to separate real stars from extended or spurious objects. 

object list (whose peak intensities were 5<r above sky back- 
ground) and to measure magnitudes both using aperture 
photometry and point spread function (psf) modeling. A 
single psf template was built using a compilation of very 
high signal-to-noise isolated stars from the different fields. 
Star/galaxy separation proceeded by applying cuts in the 
space of parameters output by the psf fitting task ALL- 
STAR. Fig. 1 shows the cuts applied to the 8m vs 7 814 
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Fig. 2. The observed composite CMD for all 32 fields in- 
cluded in our work. The magnitudes and colours are cor- 
rected for instrumental effects (aperture and CTE) and 
calibrated to the standard system. 



(panel a) and \ 2 vs Isu (panel 6) relations in the field 
uo50. From these panels we also determined cut-off mag- 
nitudes both in the bright end (due to saturation, where 
\ 2 sharply increases) and in the faint end (where dis- 
entangling stars from other sources becomes impossible). 
For the fields with too few stars for that to be done, we 
used the cut-off magnitudes from the more crowded fields, 
scaled up or down according to exposure time. 

Photometric corrections for limited aperture and 
charge-transfer effect (CTE), as prescribed by Holtzman 
et al. (1995a, b) were applied to the data. The magnitudes 
and colours were converted to the standard system, also 
using the transformations listed by those authors. Finally, 
completeness levels were obtained as a function of the 
F814TF magnitudes and used to compensate for the loss of 
stars in the faintest magnitude bins. Typically, complete- 
ness started to drop abruptly from 100% at Igi4w — 24 

In Fig. 2 we show the joint CMD for all 32 WFPC2 
fields. It contains over 1000 stars distributed within a large 
range of magnitudes and colours. The magnitudes and 
colours were obtained by aperture photometry, are in the 
standard photometric system and corrected for the several 
effects mentioned above. Notice, however, that no correc- 
tion to redenning was applied to the data. We preferred 
to incorporate redenning effects to our model predictions, 
described in the next section. 



3. Models for the Galaxy and Statistical 
Techniques 

3.1. Modeling the Galaxy 

Our main goal is to compare the observed CMD of 
HST/WFPC2 stars with theoretical CMDs based on mod- 
els for the structure of the Galaxy. The model number 
counts as a function of colour and magnitude can be 
obtained from numerical integration of the Fundamental 
Equation of Stellar Statistics. The main ingredients en- 
tering this equation are the stellar luminosity function, 
$(M) (LF), and the density distribution p(r) of stars 
for each structural component in the model. We consider 
models with 3 major components: a thin disc and a thick 
disc, both with double exponential density profiles, and 
a spheroidal or halo component, with a deprojected de 
Vaucouleurs profile. The main parameters governing these 
profiles are the horizontal and vertical exponential scales, 
ro and zq respectively, and the effective de Vaucouleurs 
radius, R e , for the halo. Also relevant are the normaliza- 
tion factors, no, for the different components, which are 
usually expressed as percentage values of the local spatial 
density of stars. 

The LF used in our modeling is the one originally 
from Wielen et al. (1983) up to My = 12. For lower 
luminosities we treat it as a power law in luminosity, 
logQ(My) oc a My with slope a considered as a free 
parameter. In our data vs model comparisons (see be- 
low), we also allow for variations in zq, ro, R e and in the 
normalization of the thick disc (non) and spheroidal com- 
ponents (nou)- For more on the models, profiles and LFs 
we refer to Reid & Majewski (1993), Majewski et al. (1996) 
and Santiago et al. (1996a). 

The numerical integration of stellar statistics was car- 
ried out for all 32 fields, taking into account their vari- 
able available magnitude ranges (set by saturation and 
detection) and completeness functions. Galactic extinction 
was also incorporated to the models using the Burstein & 
Heiles (1982) E(B-V) maps. Conversion of extinction vec- 
tors and model predictions from BV to VI was carried out 
in the same way as in Santiago et al. (1996a). The output 
from these integrals are tables with (V — I) colour and / 
magnitude counts for each field which can then be added 
together to yield a composite model CMD. 

We explored an extensive grid of models, covering large 
regions of the parameter spaces formed by the thick disc 
and halo structural parameters and LF slopes. Thick disc 
scale height and normalization values varied within the 
ranges 500 < z n < 2000 pc and 1% < n D < 8%, respec- 
tively. But notice that these two parameters were paired 
so that their anti-correlation was respected. A total of 
32 thick disc models resulted. For each one of them, we 
considered 4 choices of the horizontal thin disk c-folding 
length: r = 2500,3000,3500,4000 pc, and 2 choices for 
the LF slope at its low-luminosity end (My > 12): a = 
—0.2 and a = 0, totalling 256 thick disc models. 



4 



Kerber, Javiel & Santiago: Constraints on Thick Disc and Halo parameters 



16 



18 



20 



22 - 



24 



26 



I 1 1 1 1 I 1 1 1 1 I 1 1 1 1 I 1 1 1 1 I 

Model vs Obs. 
o Model • Obs. 



n 0D =4%; z =1200pc 
r =3500pc; flat • 





v-i 



Fig. 3. A CMD for a model realization compared 
data CMD. 
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As for the halo we used R e = 1500, 3000, 4500 pc for 
its effective radius and non = 0.1,0.125 and 0.2% for the 
local normalization. Another halo parameter varied was its 
axial ratio, with values c/a = 0.5, 0.65, 0.8, 1.0, yielding a 
total of 36 stellar halo models. 

This extensive parameters grid was defined so as to 
include the variety of parameter values often quoted in the 
literature (Norris 1999) and to accommodate most thick 
disc and halo parameters whose values are still uncertain. 

3.2. Statistical Tools and Techniques 

Fig. 3 shows a comparison of a theoretical CMD, whose 
main parameters are indicated, with the observed one. 
The simulated CMD was built from the expected model 
number counts as a function of colour and magnitude us- 
ing Monte-Carlo techniques: we randomly threw points on 
the CMD plane with probabilities proportional to the ex- 
pected number in each CMD position. The total number 
of points in the simulated CMD was a Poisson deviate of 
the total model counts. 

Our main challenge is to find an objective and efficient 
statistical tool to quantify how different or similar these 
two two-dimensional distributions of points are. A visual 
comparison reveals regions in the CMD, such as < V — 
I < 0.8 and 18 < I < 22, where there is a clear excess of 
model points. There are other regions (like < V — I < 0.8 
and 23.5 < I < 24.5, for instance) with more observed 
than model stars. The different models in the grid show a 
smooth and continuous variation in their CMDs, making 
it hard to visually select best fitting models. 



We here propose two statistics which are objective and 
simple tools for comparing two distributions of points on 
a plane. The first is simply a measure of the dispersion 
between number counts within different data bins. Let mj 
and Oi be the number of model and observed stars in the 
i th of N bins regularly spaced on the CMD plane; we then 
define the dispersion as: 



JV 



s 2 = ^(m 4 - Oif 



i=l 

The second statistics was used by Santiago & Strauss 
(1992) in a different context, to study spatial segregation 
of galaxies. It also compares the number of CMD points 
found in different bins. Let pi be the percentile position of 
\oi — mi \ within the distribution of \mi~- m,ij \ (j = 1, N rea i) 
values, where my are Poisson realizations of the expected 
model counts to; in the i th bin. We then define 



pss 



J2log(100- Pi (%)) 



We typically compute pss from a distribution of 
N rea i = 300 Poisson realizations of each model. Even 
though both statistics depend on data binning, we have 
observed that both are quite insensitive to the binning 
scheme. They are also applicable to datasets with varying 
number of points. This is in contrast with what we found 
for the W statistics (Saha 1998), whose results seem to be 
sensitive to the total number of points in both distribu- 
tions being compared. 

Clearly, the better the agreement between model and 
observed CMDs, the smaller s 2 and the larger pss will be. 
We tested both statistics by means of controlled experi- 
ments, in which we generated a Monte Carlo realization 
from a particular model (as described in the beginning 
of this section) and compared it to all the models. Fig. 4 
shows the result of a typical such experiment. It shows a 
diagnostic diagram of pss vs s 2 values obtained by com- 
paring the simulated "data" with the theoretical CMD for 
all 256 models in the thick disc grid described in the pre- 
vious section. There is a clear anti-correlation between the 
2 statistics, as expected. 

The large dot in Fig. 4 corresponds to the particular 
model used to generate the "data" points. It is clearly 
among the extreme values of both statistics, again as ex- 
pected. However, the fact that it is not the most extreme 
point in the upper left of Fig. 4 is a measure of the reso- 
lution provided by the diagnostic diagram used. 

4. Results 

4-1- Thick Disc 

In this section we apply the pss and s 2 statistics to the 
comparison between our observed CMD and the theoret- 
ical CMDs from the grid of models presented earlier. Our 
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Fig. 4. Pss vs s 2 diagram for a controlled experiment. 
Each point corresponds to a model with a particular set 
of values for the thick disc structure and LF. The large 
dot gives the position of the model from which the "data" 
CMD was created. 

main tool is the pss vs s 2 diagram. In computing pss we 
ran N rea i = 300 Poisson realizations for each model in 
order to build the distribution of rriij,j — 1,300 model 
number counts in each bin. We adopted a 14x10 binning 
covering the I vs (V-I) plane in the ranges 18 < I < 25 
and < (V - I) < 4. 

Fig. 5 shows the position of all 256 models we con- 
sidered by varying the thick disc structural parameters 
zo, noD arid tq and its stellar LF slope. There is again a 
clear correlation among the two statistics. Furthermore, 
the position within this diagram is correlated with the 
parameter values. The diagram splits into two branches, 
the upper one corresponding to models with z > 1200 
pc, which usually feature a larger than observed number 
of stars. The lower branch (small pss and s 2 values) con- 
tains models with a deficiency of stars relative to the ob- 
served sample. The best models in the upper left have 
800 < z < 1200 pc and 4 < n 0D < 8%, with little or no 
dependence on the horizontal exponential scale, r . 

These results are further confirmed by Fig. 6, where 
the pss values (assuming a decreasing LF at the faint 
end) are shown for different locations within the zq vs hod 
plane for fixed r$. The diagrams do not differ signifi- 
cantly, with the best pss values located within the range 
800 < z < 1200 pc and 4 < n 0D < 8% for all panels. 
Note that the anti-correlation between these two struc- 
tural parameters persists, a z — 800 pc with n £> ~ 8% 
being essentially as good a model as one with z — 1200 
pc and n D — 4%. 
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Fig. 5. The pss vs s 2 diagram obtained by comparing 256 
thick disc models to the real data. 
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Fig. 6. Pss maps in the z vs n D grid for 4 different 
values of the horizontal scale length r , as indicated. 



There is also a clear trend for models with decreasing 
faint end LF slope to be favoured by our observed CMD. In 
Fig. 7 we show the difference in pss values between models 
with decreasing and flat LF slopes at the faint end, all the 
other parameters being the same. The large majority of 
the 128 such differences are positive, indicating that the 
data are best described by a decreasing LF. In fact, of the 
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Fig. 7. Distribution of Apss values for pairs of thick disc 
models whose only difference is the value of a, the low- 
luminosity slope of the LF. 



35 models located closer to the upper left boundary of the 
pss vs s 2 diagram, only 1 had a flat thick disc LF slope. 
A similar result of a decreasing faint end for the LF has 
also been found for the halo and/or disk stars by Bahcall 
et al. (1994), Santiago et al. (1996a), Gould et al. (1997) 
and Mendez & Guzman (1998) using different data sets. 

The consistency of the results just shown was tested by 
adopting both coarser (70 bins) and narrower (280 bins) 
binning schemes, as well as different color ranges than the 
one quoted early in this section. There is no significant 
change in the pss vs s 2 diagram or in the position within 
it of models with different parameters. 



4.2. Halo 

In order to increase the fraction of halo stars in our ob- 
served CMD, we have considered only the 18 WFPC2 
fields with |6| > 40° and which point away from the cen- 
tre of the Galaxy. This reduced, higher latitude, sample 
contains over 200 stars. Their CMD was then compared 
to those generated by the 36 halo models discussed at the 
end of Sect. 3.1. 

The pss vs s 2 diagram for these 36 models is shown in 
Fig. 8. We again notice the anti-correlation between the 
two statistics. And again there is a correlation between 
position within this diagram and model parameters. The 
2 models above the main branch have iioh = 0.2%, c/a = 
1.0 and large values of R e (R e = 3000 or 4500 bf pc). Most 
importantly, models with c/a — 0.5 are concentrated to- 
wards the lower part of the diagram, therefore being ruled 
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Fig. 8. The pss vs s 2 diagram obtained by comparing 36 
halo models to the real data. 

out. Due to the smallness of the sample, however, it was 
not possible to single out any other region of parameter 
space favoured by the data. 

5. Summary and Conclusions 

We carried out V and I photometry for over 1000 field 
stars in the Galaxy down to faint (/ ~ 25) magnitudes 
and distributed over different directions on the sky im- 
aged with deep HST/WFPC2 exposures. The CMD for 
the data, corrected for instrumental effects and calibrated, 
was compared with theoretical CMDs expected from a grid 
of models for the structure of the Milky- Way. The main 
grid parameters are the thick disc scale height zo, its local 
density normalization, hod, the stellar halo effective ra- 
dius R e , normalization noH and axial ratio c/a, as well as 
the slope of the luminosity function at its low-luminosity 
end (M v > 12). 

We used two simple statistical tools to compare model 
and observed CMDs, in an attempt to restrict the range of 
acceptable values for the structural and LF parameters of 
the halo and thick disc. Our analysis indicates acceptable 
models in the range 800 < z a < 1200 pc, 4 < n 0D < 8% 
with these two parameters being anti-correlatcd. These 
ranges in both parameters are substantially narrower than 
those quoted in a recent review by Norris (1999). In partic- 
ular, our results do not accommodate the very thick disc 
of M > 0.3M Q stars suggested by Gyuk & Gates (1999) in 
their reinterpretation of microlensing results. The zq range 
inferred from our work is also in disagreement with earlier 
works using ground-based star counts or samples of RR 
Lyrae stars, which have yielded z > 1500 pc (Hartwick 
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1987, Rodgers 1991, Reid & Majewski 1993). On the other 
hand, as pointed out by Majewski (1993), introduction of 
metallicity gradients in our models could yield larger zq 
values. Another possibility, along similar lines, is that the 
thick disc may have a complex structure with possible 
subcomponents, such as the metal-weak thick disc (Chiba 
& Yoshii 1998, Norris 1999). Finally, our results are only 
marginally consistent with that of Robin et al. (1996), who 
find z = 760 ±50 pc and n D = 5.6 ± 1%. A better agree- 
ment is found with Buser et al. (1999), whose best model 
for the thick disc has n 0D = 5.9 ± 3% and z = 910 ± 300 
pc. 

As already established for disk stars by Santiago et 
al. (1996a), Gould et al. (1996,1997) and Mendez & Guz- 
man (1998) (from studies of HST/WFPC2 fields), we find 
the thick disc LF to be decreasing beyond My ~ 12 with 
large confidence. This is also in agreement with most stud- 
ies of globular cluster luminosity functions. 

We also attempted to restrict the structural parame- 
ters of the halo using a smaller sample with |6| > 40°. We 
were able to rule out flattened oblate halo models, with 
c/a ~ 0.5, or models with both high normalization and 
high effective radius. The smallness of the sample, how- 
ever, prevents further conclusions about halo structure. 

A way forward will be to devise more and better sta- 
tistical tools to be applied to larger stellar samples. Par- 
ticularly useful in constraining the structure of the main 
components of the Galaxy will be the inclusion of future 
kinematic data to be obtained from new orbital astromet- 
ric missions such as SIM and GAIA (Majewski 2000). 
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